/* Copyright (c) CERN
 * * Artistic License 2.0
 * */

import java.util.Random;
import java.lang.Math;

/**
 * Generalized Pareto Distribution
 * See http://en.wikipedia.org/wiki/Pareto_distribution#Generalized_Pareto_distribution
 *
 * @author Xavier Grehant
 */
public class GPD extends DistroUtils implements Distro {

/**
 * Parameter of the distribution: location.
 */
   double mu;

/**
 * Parameter of the distribution: scale.
 */
   double sigma;

/**
 * Parameter of the distribution: shape.
 */
   double xi;

/**
 * A uniform random number generator.
 */
   Random rand;

/**
 * Constructor taking the parameters of the distribution.
 * @param mu location.
 * @param sigma scale.
 * @param xi shape.
 */
   public GPD(double mu, double sigma, double xi) {
      this.mu = mu;
      this.sigma = sigma;
      this.xi = xi;
      rand = new Random();
   }

/**
 * Constructor taking the first and second moments.
 * @param xbar mean.
 * @param var variance.
 */
   public GPD(double xbar, double var) {
      mu = 0d;
      sigma = sigma(xbar, var);
      xi = xi(xbar, var);
      rand = new Random();
   }

/**
 * Constructor taking a sample set.
 * @param array the sample set.
 */
   public GPD(int[] array) {
      double xbar = mean(array);
      double var = var(array, xbar);
      mu = 0d;
      sigma = sigma(xbar, var);
      xi = xi(xbar, var);
      rand = new Random();
   }

/**
 * Constructor taking a samples set and a minimum value of interest
 * and discards values below the minimum value of interest.
 * @param array the sample set.
 * @param minvalue the minimum value of interest.
 */
   public GPD(int[] array, int minvalue) {
      int size = 0;
      for (int i = 0; i < array.length; i++)
         if (array[i] >= minvalue) size++;
      if (size == 0) construct(0d, 0d);
      else {
         int[] result = new int[size];
         int index = 0;
         for (int i = 0; i < array.length; i++) {
            if (array[i] >= minvalue) {
               result[index] = array[i];
               index++;
            }
         }
         construct(result);
      }
   }

/**
 * This method extract the parameters of the distribution from the moments
 * and define the other instance fields.
 * @return
 * @param xbar mean.
 * @param var variance.
 */
   void construct(double xbar, double var) {
      mu = 0d;
      sigma = sigma(xbar, var);
      xi = xi(xbar, var);
      rand = new Random();
   }

/**
 * This method extracts the parameters of the distribution from a sample set
 * and define the other instance fields.
 * @return
 * @param xbar mean.
 * @param var variance.
 */
   void construct(int[] array) {
      double xbar = mean(array);
      double var = var(array, xbar);
      construct(xbar, var);
   }

/**
 * This factory method returns a biased GPD estimation to fit samples over minvalue.
 * @return a GPD.
 * @param values sample set.
 * @param minvalue minimum value of interest.
 * @param refinements number of times an estimation is biased to better fit over
 *        the minimum value of interest.
 */
   static GPD biased(int[] values, int minvalue, int refinements) {
      GPD gpd = new GPD(values);
      for (int refined = 0; refined < refinements; refined++) {
         values = gpd.bias(values, minvalue);
         gpd = new GPD(values);
      }
      return gpd;
   }

/**
 * Alpha
 * Used in the random number generator.
 * @param xbar mean.
 * @param var variance.
 */ 
   double alpha(double xbar, double var) {
      return .5d * xbar * (Math.pow(xbar, 2) / Math.pow(var, 2) + 1d);
   }

/**
 * Kappa
 * Used in the random number generator.
 * @param xbar mean.
 * @param var variance.
 */
   double kappa(double xbar, double var) {
      return .5d * (Math.pow(xbar, 2) / Math.pow(var, 2) - 1d);
   }

/**
 * Sigma
 * Used in the random number generator.
 * @param xbar mean.
 * @param var variance.
 */
   public double sigma(double xbar, double var) {
      return alpha(xbar, var);
   }

/**
 * Xi
 * Used in the random number generator.
 * @param xbar mean.
 * @param var variance.
 */
   public double xi(double xbar, double var) {
      return -kappa(xbar, var);
   }

/**
 * Random number generator according to the distribution.
 * See Handbook of statistical distributions with applications
 * Krishnamoorthy, K.
 * Chapman and Hall/CRC Press, 2006
 * @return random number according to the distribution.
 */
   public double nextDouble() {
      return mu + ((sigma * (Math.pow(rand.nextDouble(), -xi) - 1d)) / xi);
   }

/**
 * Cumulative density function of the distribution.
 * @return P(X <= 'limit')
 * @param the point at which the cdf is calculated.
 */
   public double cdf(int limit) {
      return 1d - Math.pow(1d + (xi * (double) (limit - mu) / sigma), - 1d / xi);
   }
}
